Molecular simulation and theory of a liquid crystalline disclination core 
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Molecular simulations of a nematic liquid crystal confined in cylinder geometry with homeotropic 
anchoring have been carried out. The core structure of a disclination line defect of strength +1 has 
been examined, and comparison made with various theoretical treatments, which are presented in 
a unified way. It is found that excellent fits to the cylindrically-symmetrized order tensor profiles 
may be obtained with appropriate parameter choices; notwithstanding this, on the timescales of the 
simulation, the cylindrical symmetry of the core is broken and two defects of strength +1/2 may be 
resolved. 

PACS numbers: PACS: 61.30.Cz,61.30. Jf,61.20. Ja,07.05.Tp 
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I. INTRODUCTION 



Liquid crystals are at the heart of a range of techno- 
logical devices, which rely on the ability to manipulate 
the direction of preferred molecular alignment, the di- 
rector n, with electric and magnetic fields, and through 
coupling to surfaces. The practical performance of such 
devices, and the theoretical description of director distor- 
tions, both rely on smooth variation of the director with 
position in space; this is well accounted for by the Frank 
free energy g], | which is an expansion in squared gradi- 
ents of the director field n(r), parametrized by the splay 
{Kn), twist (K22) and bend (-K33) elastic constants. The 
director n is the principal axis of the local second-rank 
ordering tensor Q, which characterizes the nematic state; 
such a description is intrinsically uniaxial, and neglects 
variation of the degree of ordering with position, of the 
kind which occurs near bounding surfaces and around 
topological defects. Such defects are treated as singu- 
larities in the director field; for a better description, it is 
necessary to replace the director field by one which allows 
some variation of the relevant order tensor components 
over relatively short length scales. The relevant region 
near the defect is called the core. 

The most common type of defect in the nematic phase 
is the disclination line defect, characterized by an integer 
or half-integer index defined by the number of turns of 
the director field associated with taking a circuit about 
the line. The phenomenological approach to the inves- 
tigation of line defects was described by Schopohl and 
Sluckin ||, . They worked with the full order-parameter 
tensor Q a p in the framework of Landau-de Gennes the- 
ory, and applied this to the structure of the ±^ discli- 
nation core. The full order parameter Q a p was used a) 
to avoid divergent terms in the elastic energy; b) to take 
into account possible biaxiality of the defect core. Then, 
the general idea was implemented for the particular prob- 
lem of the core structure of the ±1 strength defect: Bis- 
cari and Virga Q and Mottram and Hogan |6| solved 
the equations for the order tensor and obtained order 
parameter and biaxiality profiles. They used truncated 
expansions of the free energy density, which helped them 
to obtain analytical solutions to the problem. Finally, 



Sigillo et al considered the +1 disclination problem in 
the spirit of Maier-Saupe mean-field theory. 

In this paper we use the approach proposed by 
Schopohl and Sluckin to obtain the order parameter pro- 
files near the disclination line. We also recapitulate the 
results of the other theories in order to discuss the ad- 
vantages and disadvantages of the proposed models. We 
compare these predictions with the results of computer 
simulation using a simple molecular model. Computer 
simulation is a well-established method of relating bulk 
elastic coefficients j|, [| [t(| [tl], [l2[ O, |1 and more re- 
cently surface anchoring strengths fig , 16 1, to molecular 
structure and interaction parameters. An early study of 
disclination line defects |lj involved Monte Carlo simu- 
lations of rod-like molecules (hard spherocylinders) con- 
fined in cylindrical geometry with boundary conditions 
chosen to stabilize the chosen director field far from the 
defect. This work concentrated on the disclination line 
defect of strength — |, in cylinders of radius 2-3 times 
the molecular length. For short rods, smooth variation 
of the order parameters with position was observed, and 
the defect core retained axial symmetry. For longer rods, 
having much larger bend elastic constants, this symmetry 
was broken, and various microscopically 'escaped' struc- 
tures were seen. There was also evidence of metastability, 
and non-convergence of structures from different starting 
configurations, which the authors attempted to resolve 
by free energy calculations. More recently, a study of 
two-dimensional models has been carried out |l8| . These 
studies support the view that the disclination core is one 
or two molecular lengths across. 

The purpose of this paper is to present simulation re- 
sults for a model of hard particles in a cylindrical pore, 
following closely the approach of Ref. |l7j , and compare 
with the theories just mentioned. We study the simplest 
example of biaxial molecular arrangement, namely the 
disclination defect of strength +1 corresponding to a uni- 
form, cylindrically symmetric, splay deformation of the 
director far from the core. This is imposed by choosing 
homeotropic anchoring conditions at the cylinder pore 
walls. 

We take care to check the equilibration of our simula- 
tions, and note that quite extensive runs are required to 
ensure this. We also pay particular attention to studying 
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the effects of pore size, with radius varying from 2 to 5 
times the molecular length, so as to eliminate the effects 
of the walls on the defect core. More precisely, our aim 
is to restrict these effects to those which arise from the 
planar radial far-field structure, rather than from density 
and order parameter variations in the immediate vicin- 
ity of the walls. It is well-known that the pore radius 
plays a critical role in determining the stability of dif- 
ferent nematic structures in this geometry. An analysis 
of the elastic free energy shows that the 'escaped radial' 
structure, in which the director bends over to become 
parallel with the symmetry axis in the core region, is the 
stable structure for large pore radius [[19], ^fj. In analyz- 
ing experiments on nematic liq uid crystals in cylindrical 
pores, Crawford et al. ^2| show that at small pore 
radii, the escaped structure (which will also, in general, 
have point defects along the disclination line) is not the 
most stable; they discuss the planar radial and planar 
polar configurations, finding the latter to be stable for 
the parameter values that they survey. 

A unified presentation of the different theories is given 
in Sec. 55 the model, and simulation techniques, are set 
out in Sec. |l|. The observed structures are described in 
Sec. |y], along with theoretical fits to the order parame- 
ter profiles. A discussion of the results, and the validity 
of the theories, together with some concluding remarks, 
appear in Sec. [v|. 



II. PHENOMENOLOGICAL MODELS 



In the framework of the continuum theory, the system 
can be described by the Landau-de Gennes free energy 
density @: 



F = k | VQ| + a (Q) 



(1) 



where a (Q) is a function of the invariants of Q, the 
symmetric order tensor. This is the suitably-normalized 
traceless part of the tensor M of the second moments of 
the molecular orientational distribution function / 



M = J /(u) u®ud!!, 



(2) 



where I is the unit tensor and the integration is over the 
unit sphere. 

We exclude the escape of the director in z direction, 
and our boundary conditions also prevent spiral config- 
urations, so the eigenvectors of Q coincide with the unit 
vectors of the cylindrical geometry, namely e p , eg and 
e z , respectively, in the radial, tangential, and axial di- 
rections. Since the eigenvalues of Q are not independent, 
the tensor may be expressed in terms of two independent 



parameters: the order parameter S and biaxiality a 
Q = Q PP e p (g) e p + Qee e g ®e g + Q zz e z <g> e z 
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(3) 



Minimizing the full free energy functional ([!]) with re- 
spect to the S, a and taking into account that S = S (p), 
a = a (p) , where p is the distance from the cylinder axis, 
we obtain the following Euler-Lagrange equations: 



[pS')' 3 P 
(pa')' + p 



- 1 (S-a)-3g 1 (S,a) = 0, 
^(S-a)-g 2 (S,a) = 0. 



(4) 



Here g\ (S, a) = ^K~ 1 da/dS, g 2 (S,a) = ^n~ 1 da/da, 
and the prime denotes differentiation with respect to p. 

The system of equations (Q) should be solved with cor- 
responding boundary conditions. The boundary condi- 
tions at p — can be obtained if we look for the solutions 
close to the centre of the disclination line. Indeed, in the 
limit p — > 0, we can neglect the functions 31,32 in equa- 
tions (|4|). Then, seeking the solution as an expansion in 
powers of p, one can obtain that in the region close to 
the centre of the disclination line the solutions are 



S = S + 37p 2 , a = S - ■yp 2 



(5) 



where So and 7 are constants. Therefore, at p = 0, we 
have 



p=0 



. 



n 



p=0 



. 



(G) 



We also assume that, far away from the disclination core, 
we have uniaxial nematic, and the surface at p = R pro- 
vides the order parameter S s : 



S\ 



p=R 



5', 



\p=R 



. 



(7) 



Some general properties of the equations (^), regard- 
less of the explicit form of the functions 31, 32, can help 
us to fit the simulation results. The simulations provide 
the value of the order parameter on the disclination line 
So which we can also derive analytically. Indeed, the 



solutions (H) imply that S + 3a 



p^O 



const. The con- 



dition that the energy is bounded requires also So = cto- 
Therefore, from (Q) we obtain the implicit equation for 
the order parameter value on the disclination line: 

31 (So, S ) +92 (So, S ) = 0. (8) 
Now we consider the particular models. 

A. Full free energy expansion 

A widely used formula for a is Landau-de Gennes': 

(J lg = aTrQ 2 - 6TrQ 3 + c [TrQ 2 ] 2 , (9) 
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where a is assumed to depend linearly on the tempera- 
ture, whereas positive constants b, c are considered tem- 
perature independent. For this free energy, the uniaxial 
nematic state is stable when b 2 > 2Aca with the degree 
of orientational ordering 



This potential provides the following functions gi,g2- 



S b 




(10) 



Taking into account the parametrization (|3j) the suitably 
rescaled potential (ph may therefore be rewritten as 



= { ^ ( \s 2 - \s (S b + S U ) + S b S u 



3S u S b + 6(S U + S b )S + 3S 2 + -a 2 



(11) 



so that the turning points for the uniaxial nematic phase 
with a — occur at S = S u , S b . 
The functions g\ , gi then read 

\ 2 9l = s u S b S + (S u + S b ) {3a 2 - S 2 } + S {3a 2 + S 2 } , 
A 2 <?2 = 3S u Sb<x + 6 (S u + Sb) Sa + 3a {3a 2 + S* 2 } , 

where A 2 = jK/k is the characteristic length. 

The equation (||) provides us with the values of the 
order parameter and biaxiality on the disclination line: 



So = a 



1 



:S U 



(12) 



B. Biscari and Virga approach 



In order to obtain analytical expressions for the order 
parameter profile, Biscari and Virga JE) used a quadratic 
approximation to the full free energy a, expressed in 
terms of S and a: 

<J BV = k { V (S - S b ) 2 + a 2 } (13) 
with the following expressions for the functions gi,g2- 

gi = (S - S b ) , .92 = ^« • 

The order parameter at the center of the disclination line 
then reads 



So 



>1 



1+7] 



- Si, 



(14) 



C. Mottram and Hogan approach 

Another model, considered by Mottram and Hogan || 
uses a quartic potential in S but retains a quadratic po- 
tential in a 

omh = k Qs 2 - ^S (S b + S u ) + SbSu^j + a 2 

(15i 



.91 = -j^S (S - S u ) (S - S b ) 



92 = 



and the following value of the order parameter on the 
disclination line: 



S = i ls u + S b ± s] (S u - Sbf - 47T 1 ] (16) 



D. Maier-Saupe approach 

Another mean-field approach to the description of the 
disclination line was considered by Sigillo et al . They 
applied Maier-Saupe theory and obtained expressions for 
the functions 51,32- The Maier-Saupe theory has an in- 
termolecular potential strength U as the only adjustable 
parameter, determining the correct value of the bulk or- 
der parameter Sb (U). The value of the order parameter 
So on the disclination line is then fixed and determined 
by the value of Sb- 



III. MOLECULAR MODEL AND SIMULATION 
METHODS 

The molecules in this study are modelled as hard ellip- 
soids of revolution of elongation e = A/B = 15, where A 
is the length of the major axis and B the length of the two 
equal minor axes. The phase diagram and properties of 
this family of models are well studied @ || ||, ||, |§. 
Units of length are chosen such that AB 2 = 1, making 
the molecular volume equal to that of a sphere of unit 
diameter. It is useful to express the density as a frac- 
tion of the close-packed density for perfectly aligned hard 
ellipsoids, assuming an affinely-transformed face-centred 
cubic or hexagonal close packed lattice; in reduced units 
g* p = Q C pAB 2 — \/2. Henceforth the asterisks denoting 
reduced quantities will be omitted. 

The molecules were confined within a cylinder of ra- 
dius i?, and height H, with periodic boundary condi- 
tions applied in the z direction (the symmetry axis). 
The cylindrical confining walls are defined by the con- 
dition that they cannot be penetrated by the centres of 
the ellipsoidal molecules. Packing considerations gener- 
ate homeotropic ordering at the surface, i.e. the molecules 
prefer to orient normal to the interface, without the need 
to apply an explicit ordering field. The properties of 
such surfaces in planar geometry were investigated pre- 
viously |]l6| . For these particles, the isotropic-nematic 
phase transition occurs at quite low density, g/ g cp « 0.2. 
Temperature is not a significant thermodynamic quantity 
in this model; we set fceT = 1 throughout. 

Monte Carlo simulations were carried out for sys- 
tems of the following cylinder radii: R/A = 
2.08,2.67,4.00,5.33. The corresponding numbers of el- 
lipsoids were N = 3500, 6000, 13000, 22000. The height 
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H/A was in the range 2.6-2.8 in all cases. All the sim- 
ulations were conducted at the same state point used in 
the earlier study corresponding to a bulk pressure 
P = 2.0 in the above reduced units; this is significantly 
higher than the isotropic-nematic coexistence pressure 
P ss 1.49 [^tJ. No external fields were applied, and con- 
ventional Monte Carlo moves were employed, with trans- 
lational and rotational displacements chosen to give a 
reasonable acceptance rate. 

The simulation results were analyzed to give a density 
profile, averaged over cylindrical shells of width 0.125B- 
0.25-B, and an order tensor profile obtained by averaging 

Q a p(pi) = ^^UjaUjp - 5<wj^ at,f3 = p,6,z 

where there are molecules present in bin i; 8 a p is the 
Kronecker delta. The axis system is resolved as before 
into radial (p), tangential (#), and axial (z) components, 
for the purposes of accumulating these functions. There- 
fore, the tensor components Q a p have been averaged over 
global rotations of the system (positions and orientations 
together) about the symmetry axis, as well as global 
translations in the z direction. Diagonalizing this ten- 
sor, for each bin, gives three eigenvalues and three corre- 
sponding eigenvectors. This procedure allows us to test 
for (a) the planar radial structure, for which the eigen- 
vectors lie along the p, 9, z coordinate axes, independent 
of p, and the eigenvalues are Q pp (p), Qee{p), Qzz(p)); 
and (b) the escaped radial structure, for which one of 
the eigenvectors points along the 8 direction for all val- 
ues of p, and the other two lie in the p, z plane, changing 
their orientation as p varies. 

The above tensor is not sensitive to any breaking of 
axial symmetry: instead, it gives a cylindrical average. 
To give a full representation of the positional variation 
of the orientational order, we must retain the full spatial 
dependence of Q. We find it convenient to calculate 

S a p(r) = — u ja Ujf} a,f3 = x,y,z (17) 

where the sum is conducted over the n-jz molecules found 
to be in a neighbourhood 1Z of the chosen point r. The 
eigenvectors of this S tensor are the same as those of 
the corresponding Q; the eigenvalues are linearly related 
and are non-negative, so the tensor may be visually repre- 
sented as a spheroid, whose principal axis directions and 
lengths are given by the eigenvectors and corresponding 
eigenvalues p8|. Then, a well-ordered uniaxial region ap- 
pears as an elongated, prolate, shape; the uniaxial core 
of a disclination defect of strength +1 would appear as 
an oblate spheroid with its symmetry axis along the axis 
of the cylinder; disordered regions would correspond to a 
spherical shape, and so on. 
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FIG. 1: Time evolution of system started from perfect pla- 
nar radial configuration, in cylinder of radius R/A = 2.67. 
We show order tensor eigenvalue profiles averaged over 10 5 
sweeps, taken at the indicated times (units of xlO 6 sweeps). 
The top eigenvalue is Q pp , the middle one Qee and the bot- 
tom one Q zz in each case. In the upper panel we plot the 
biaxiality parameter a — \{Qee —Qzz)- Note that the biaxial 
core develops quite rapidly, and that equilibration is essen- 
tially complete between 0.4 and 0.8 xlO 6 sweeps. The final 
equilibrium configuration is identical with that obtained from 
a disordered starting configuration. The shaded ellipse indi- 
cates the outermost, homeotropically anchored, layer at the 
cylinder wall. 



IV. SIMULATION RESULTS 

Most results were obtained from starting configura- 
tions in which the ellipsoids were positioned randomly, 
avoiding overlaps, but perfectly aligned in the planar ra- 
dial configuration. Typical equilibration times for these 
systems were around 7 x 10 5 MC sweeps (one sweep is 
one attempted move per particle), and this was mon- 
itored through the orientation profiles defined above. 
The development of a biaxial core in the system with 
R/A = 2.67 is shown in Fig. [jj the other systems evolved 
in a similar way. Following equilibration, production runs 
of approximately 5 x 10 5 MC sweeps were undertaken. 

Checks for convergence were carried out, starting from 
configurations in which molecular orientations as well as 
positions were disordered. The time evolution of the 
R/A = 2.67 system is shown in Fig. 0. Note that the or- 
dering remains uniaxial at the boundary, and propagates 
in towards the centre; as the ordering reaches the centre, 
the biaxial disclination core develops. An equilibration 
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FIG. 2: Time evolution of system started from orientation- 
ally disordered configuration, in cylinder of radius R/A = 
2.67. Notation as for Fig. [|. Note that the ordering remains 
uniaxial at the boundary, and propagates in towards the cen- 
tre; as the ordering reaches the centre, the biaxial disclination 
core develops. The final equilibrium configuration is identical 
with that obtained from a planar radial starting configura- 
tion. 



period of 10 6 sweeps was needed to reach the equilib- 
rium planar radial structure, and a subsequent produc- 
tion run of 2 x 10 5 sweeps yielded identical results to 
those obtained from the planar radial starting point. For 
the larger cylinder radii, similar behaviour was observed, 
but the timescale for propagation of the order from the 
boundary to the centre was correspondingly longer. In 
each case, a run was conducted which was of sufficient 
length to confirm that the correct structure was becom- 
ing established. 

In no case was the escaped radial structure formed dur- 
ing these runs. We have carried out some preliminary 
tests, for the pores of smaller radius, R/A — 2.08, 2.67, in 
which the escaped radial structure was stabilized by the 
application of a uniform orienting field favouring align- 
ment along the cylinder axis. Following removal of this 
field, the planar structure was seen to be recovered on 
a simulation timescale of 5 x 10 5 sweeps. Thus, the pla- 
nar structure is thermodynamically stable for these cases, 
and it appears to be at least metastable on the timescales 
of our simulations for the larger pores. 

In Fig. ID the order parameter profiles are presented 
for all the different cylinder sizes. The biaxiality profile 
a(p) — \(Qee ~ Qzz) indicates the extent of the core 
region. For the smallest two cylinders, it is clear that 



FIG. 3: Equilibrium order tensor eigenvalue profiles for cylin- 
ders of indicated radii. For the smallest two cylinders, the 
walls clearly have an effect on the defect core structure; the 
results for the two largest cylinders are essentially indistin- 
guishable. Notation as for Fig. |l| 



the walls act to deform the disclination core. For R/A = 
2.08, the core region is enlarged, while for R/A = 2.67 it 
is compressed. The results for the two largest cylinders 
are almost indistinguishable. 

In Fig. |] we show the order tensor variation in the 
xy plane (no variation with z-coordinate was observed) 
for the core region in the R/A = 4.00 case. The S- 
tensor spheroid of eqn (|l7|) is plotted at a number of 
points (x s ,y s ) which lie on circles centred on the cylin- 
der axis, separated by 0.25A, with the neighborhood 
1Z of each point defined to include all molecu les whose 
centres (x, y, z) satisfy \J {x — x s ) 2 + (y — y s ) 2 < 0.25A 
The tensors are averaged over a run of 10 5 sweeps. The 
figure indicates that the core is actually best described 
as a pair of disclinations each of strength +1/2, sym- 
metrically arranged at p/A «1-1.25 from the cylinder 
axis. Very similar results are seen for all the three largest 
pores, R/A — 2.67, 4.00, 5.33, while for the smallest pore, 
R/A — 2.08, the defects are slightly further out, at 
p/A ps 1.5. There is a small region of almost unper- 
turbed uniaxial nematic liquid crystal around p ~ 0, with 
the director perpendicular to the cylinder axis. We note 
that similar defect pairs are seen in the two-dimensional 
simulations reported in Ref. Jl8[ . 

The order tensor profiles reported in Fig. ^ are prop- 
erly regarded as axial averages of structures which are 
not themselves axially symmetric on the timescales of 



FIG. 4: Distribution of order tensor in the plane perpen- 
dicular to the cylinder axis, in the core region, for the case 
R/A = 4.00. Concentric circles at half-integer intervals of ra- 
dius p/A are plotted as a guide. The spheroid sizes are not 
physically significant. The figure is oriented so that the +1/2 
defects lie above and below the cylinder axis. 



the simulation. This is why the lowest eigenvalue in the 
profiles of Fig. || adopts similar values inside and out- 
side the core: it corresponds to the eigenvector along the 
cylinder axis, and is unaffected by the axial averaging. 
We discuss in section [v| the validity of carrying out such 
an axial average, simply noting here that it is the most 
straightforward way of comparing our simulation results 
with theories which assume cylindrical symmetry. 

We have fitted the results for R/A — 4.00 using the 
phenomenological theories of section O. In each case, the 
boundary- value problem of eqns (Q) , (pj) , (0) , was solved 

||, the fitting 



using the relaxation method 29 




curves and the original simulation data are shown. Note, 
that to fit the simulation data we first determind the 
value of the parameter rj or S u using equations (|T^), (14), 
( |l6| ) . Then we performed suitable rescaling, changing the 
factor A. 



V. DISCUSSION 

One can see in Figs ^, ||, that the slope of the fit- 
ting curves in general reflects the structure of the core: 
the center of the core is strongly biaxial, extending over 
a few units of length, as shown by the splitting of the 
eigenvalues Qee, Qzz, and hence the non-zero biaxiality 
parameter a(p). 

At the same time, the difference between the descrip- 
tions is also evident. Biscari and Virga's theory gives 
an incorrect shape of S (p) for small values of S (and 
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FIG. 5: Fits of simulation results (solid lines) to theoretical 
predictions discussed in the text. The best fit in each case was 
obtained for the following parameters: Mottram, Hogan |]] 
(dotted lines): eqn ((TH)), with n/X 2 = 1.58, A" 2 = F R = 6.5, 
S b = 0.88 S s = 0.82, Su = 0.01. Biscari, Virga [§] (dashed 
lines): eqn @, with rj/X 2 = 1, A -2 = 3.4, R = 6, S b = 0.88, 
S s = 0.82. Sigillo, Greco, Marrucci |R (long dashed lines): 
U = 13, S s = 0.82, R = 7. Notation as for Fig. 0. 



hence, small values of p here). This is fairly predictable, 
since a quadratic expansion of the free energy density 
F ~ i] (S — Sb) 2 was used, which is valid only for small 
deviations of the order parameter from the bulk value Sb- 
The more sophisticated form of the free energy used by 
Mottram and Hogan (up to fourth order in S) corrects 
the slope of the curve for the order parameter S = Q pp - 
However, the biaxial part, Qgg — Q zz , still has only qual- 
itative agreement with the simulation results, which is 
probably because of the quadratic approximation to the 
biaxial part of the free energy. This is particularly ap- 
parent in the a(p) profile, in which the biaxiality extends 
far beyond the core radius which would be deduced from 
Qpp- 

In spite of having only one adjustable parameter, the 
Maier-Saupe theory of Sigillo et al. gives a very realistic 
description of the disclination core structure. However, 
it predicts slightly lower core biaxiality than obtained in 
the simulations. 

The most accurate fitting we obtained used the full ex- 
pansion of the free energy in powers of the order param- 
eter tensor (^]), as shown in Fig. ^. This form, originally 
used by Schopohl and Sluckin, manages to reproduce the 
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FIG. 6: Fit of simulation results (solid lines) to the Landau- 
de Gennes free energy expansion ( [ill ) using the approach of 
Schopohl and Sluckin JiJ ^| (dashed lines). The best fit in 
each case was obtained for A~ 2 = 1.42, R = 6, Sb = 0.88, 
S s = 0.82, S u = -0.4. Notation as for Fig. [l[ 

overall extent of the biaxial region quite well, while fitting 
the limiting behaviour of both a{p) and S(p) as p — > 0, 
and the magnitude of S in the bulk. Therefore, we can 
conclude, that all the terms in the free energy expansion 
(or, at least up to the fourth order in the order param- 
eter tensor) should be taken into account if one want to 
make the correct quantitative description of the disclina- 
tion core structure. 

Fig [| shows that the core structure seen on the simu- 
lation timescale is not, in fact, cylindrically symmetrical, 
and this raises questions regarding the validity of aver- 
aging the order tensor over rotations about the cylinder 
axis. Such averaging may arise in a natural way. We 



see a small amount of rotation of the defect pair around 
the axis, of the order of 10-15 degrees, as the core struc- 
ture evolves in time, during the 5 x 10 5 -sweep simula- 
tion runs conducted here. If we roughly equate this to 
a real-time period of the order of nanoseconds, it seems 
possible that significant rotation will occur on the exper- 
imental timescale. The situation is further complicated 
when one considers correlations along the z-direction; the 
periodic box employed here is quite small, but use of a 
much longer cylindrical pore might reveal some twisting 
of the separate +1/2 defect lines about the cylinder axis. 
Both effects would result in the kind of cylindrical aver- 
aging which we have carried out, but we can currently 
only speculate on this. 

A second question concerns the distance between the 
+1/2 defect lines. We do not observe separation of the 
defect pair, either as the simulations proceed in time, or 
as we study larger pores. Nonetheless, it is quite possible 
that the pore walls are confining the defect pair to the 
vicinity of the axis, and that they might separate if a 
much larger system were simulated. We note that, if they 
were to completely separate and approach the pore walls, 
the result would be the planar polar structure discussed 
elsewhere pi], [22] ]. It would be interesting to calculate 
free energies as a function of the defect separation to 
investigate this. 

Finally, we may expect the escaped radial structure 
to become increasingly favoured as the pore radius in- 
creases; it is clearly not the most stable structure for 
the smallest pores studied here. Once more, free energy 
calculations are needed to make a proper comparison. 
Further work on these aspects is in progress. 
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